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Abstract 



. The coupled dynamics of the isovector and isoscalar giant quadrupole resonances and 

I low lying modes (including scissors) are studied with the help of the Wigner Function 

■ Moments (WFM) method generalized to take into account pair correlations. Equations 

of motion for collective variables are derived on the basis of the Time Dependent Hartree- 

Fock-Bogoliubov (TDHFB) equations in the harmonic oscillator model with quadrupole- 
O ; 

^ I quadrupole (QQ) residual interaction and a Gaussian pairing force. Special care is taken 



of the continuity equation. 



1 Introduction 



An exhaustive analysis of the dynamics of the scissors mode and the isovector giant quadrupole 
00 ' 

^ I resonance in a harmonic oscillator model with QQ residual interaction has been performed in 

[T]. The WFM method was applied to derive the dynamical equations for angular momentum 
and quadrupole moment. Analytical expressions for energies, B{M1)- and i?(i?2)-values, sum 
rules and flow patterns were found for arbitrary values of the deformation parameter. These 
calculations were performed without pair correlations. However, it is well known [2] that pair- 
ing is very important for the correct description of the scissors mode. A flrst attempt to include 
pairing into the WFM method was made in [3j, where the description of qualitative and quan- 
titative characteristics of the scissors mode was drastically improved. However, the variation of 
the gap during vibrations was neglected there, resulting in a violation of the continuity equation 
and in the appearance of an instability in the isoscalar channel. In the present work we suggest 
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a generalization of the WFM method which takes into account pair correlations conserving the 
continuity equation. 



2 Phase space moments of TDHFB equations 

The time dependent HFB equations in matrix formulation are [5', '5] 

tnn = [n,n] (i) 

with 

The normal density matrix p and Hamiltonian h are hermitian whereas the abnormal density 
k and the pairing gap A are skew symmetric: = —k*, A"!" = —A*. 
The detailed form of the TDHFB equations is 

ihp = hp ~ ph — Ak^ + kA^, —ihk = —hk — kh* + A — Ap* — pA, 
-ihp =h*p* -p*h* -A^k + k^A, -ihfz =h*k^ + k^h-A^ + A^p + p*Al (3) 

We will work with the Wigner transformation [S] of these equations. The relevant mathematical 
details can be found in [3]. From now on, we will not specify the spin and isospin indices in 
order to make the formulae more transparent. The isospin indices will be re- introduced at the 
end. In addition, we will not write out the coordinate dependence (r, p) of all functions. The 
Wigner transform of ([3]) can be written as 

ihf = ih{h,f}-AK* + KA*-^-^{A,K*} + ^-^{K,A*} 
-j[{{n,A*}}-{{A,K*}}] + ..., 

-ih'f = ih{hJ}-A*K + K*A-^-^{A*,K} + j{K*,A} 

A*}} -{{A, «:*}}] + ..., 
—ihk = —hn — nh — —{/;,, k} — —{/€, /i} 

+A-A/-/A-|{/,A}-|{A,/} 



—ihk* 



+y [{{/^, '^D + {{'^, h}} + {{A, /}} + {{/, A}}] + 

K*h + JiK* + h} + ^ih, K*} 

-A* + /A* + A7 + y {/, A*} + y {A*, /} 

- vKi'^*' M} + {{/^, '^*}} + {{/, A*}} + {{A*, /}}] + 



(4) 



where the functions h 



/, A, and k are the Wigner transforms of h, p, A, and k, respectively, 



/(r, p) = /(r, — p), {/, g} is the Poisson bracket of the functions /(r, p) and g{T, p) and {{/, g}} 
is their double Poisson bracket; the dots stand for terms proportional to higher powers of h. 

In order to study collective modes described by these equations, we apply the method of 
Wigner function moments (or phase space moments). The idea of the method is based on 
the virial theorems by Chandrasekhar and Lebovitz [6J; its detailed formulation can be found 
in [71 IB]- For the investigation of the quadrupole collective motion with = 1+ in axially 
symmetric nuclei, it is necessary to calculate moments of Eqs. @ with the weight functions 



This procedure means that we refrain from seeking the whole density matrix and restrict our- 
selves to the knowledge of only several moments. Nevertheless, this information turns out to 
be sufficient for a satisfactory description of various collective modes with quantum numbers 
K"^ = 1+, as it was shown in our previous publications [H [3, [8]. In the case without pairing, 
this restricted information can be extracted from the TDHF equations and becomes exact only 
for the harmonic oscillator with multipole-multipole residual interactions. For more realistic 
models it becomes approximate even without pairing. The TDHFB equations (jlj) are consider- 
ably more complicated than the TDHF ones, so additional approximations are necessary even 
for the simple model considered here, as will be discussed below. 

Integration of Eqs. (jlj) (including the terms of higher orders in h) over the phase space with 
the weight W, where W is any one of the weight functions listed in ([5]), yields the following set 
of equations: 



W 



xz, p^pz, zp^ + xp^ = L, and zp^ - xp^ = /, 



(5) 
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i{{W,A*}}K-{{W,A}}K*) 




d{p,r)Wf 
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{{{W,A*}}^-{{W,A}}^*) 




d{p, r)WK 





+ '-^{W, A}(/ - /) - j[{{W, {h + h)}}^ + {{W, A}}(/ + /)] 



+j{W, A*}{f - /) + ^[{{W, {h + h)}}^* + {{W, A*}}(/ + /)]] , (6) 



where Jd{p,r) = 2{2'n-h)^^ Jd^p Jd^r. It is necessary to note an essential point: there are no 
terms with higher powers of h in these equations. The infinite number of terms proportional 
to h"" with n > 2 have disappeared after integration. This fact does not mean that higher 
powers of h are not necessary for the exact solution of the problem. As it will be shown below, 
the equations contain terms which are coupled to dynamical equations for higher-order 
moments, which include, naturally, higher powers of h. 

2.1 Continuity equation 

The moments of the first two equations with the weight xz are of special interest, because in 
this case equations are integrated over momentum space without any weight. It is known that 
the integration of the Vlasov (or Boltzmann) equation over p gives the continuity equation [5] ; 
the same is true for the Wigner function equation. The inclusion of pair correlations must not 
destroy this property of the Wigner function, so one can expect that integration over p of the 
first (and second) equation in (j4]) will produce the continuity equation. Though this is known 
from general arguments, let us repeat it in detail. 

First of all we integrate over p the first part of the first equation in (jl]): 




(terms containing A and n). 



(7) 
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In the case of a velocity-independent potential, the first integral on the right-hand side is equal 
to zero. By definition Jd^pf{r,p) = n(r) and Jd^ppif{r,p) = mn(r)-Ui(r), where n{r) is the 
density of particles and Mj(r) is the i-th component of their mean velocity. Using dh/dpi = Pi/m, 
one obtains 

n + div{nu) = terms containing A and k. (8) 

So, we will recover the continuity equation if the integral over p of all the terms of the first 
equation in (jl]) containing A and k gives zero. 

Let us integrate over p the term kA* — k*A. It is known that the pairing gap and the 
abnormal density are connected by the integral relation (see, e.g., ref. [5]): 

d^p' 



A(r, P) = - / (^^(IP - P'l)'^(r, P')- (9) 
Using this relation, one finds 

J d^pi^A* - K*A) = Jd^pJ d^p'[K{r, p)K*(r, p') - K*{r, p)K(r, p )]^^(|p - p'|) = 0. (10) 

The equality becomes obvious after changing the variables p ^ p' in the second (or first) part 
of this formula. 

Finally, let us integrate over p the term with Poisson brackets: 



d'pjd'p'j: 

1=1 



3 ^ dK*(v,p)dv{\v-v'\) g/t*(r,p) n. d<^.v') 

dxi dpi ' dpi dxi 

(9K(r,p)(9t;(|p-p'|) , <9K(r,p) 9K*(r,p') 
+ — ^ ^ (r,p) ^ t;(|p-p I)- 



dxi dpi ' dpi dxi 



d'pjd'p'j: 

1=1 



K[r,p) + K (r,p) 



oxj ox. 



, dK{r,p) dK*{r,p') 



^%^ = 0. (11) 
dpi 



dxi dxi 
The last equality becomes obvious after changing the variables p *^ p' and using the relation 
dv{\p — p'\)ldpi = —dv{\p — p'D/dp'j^. In a similar way one can show that integration over p 
of any term of higher order in h will give zero. 

So, as expected, one can conclude that pairing does not spoil the continuity equation, which 
is contained in the TDHFB equations. As a result, the A-dependent terms in the first and 
second equations of disappear, when the weight W does not depend on p, for example, 
W = xz (see below). 



2.2 Linearization 



It is convenient to rewrite the equations ([6]) in terms of h± = h ± h, f± = f ± f, Aj- = 
A± A*, K,± = K±K*. These equations are strongly nonhnear. Having in mind small amplitude 
oscillations, we will linearize: f± = f±+Sf±, k± = k,^+6k,±, A± = A^+6A±. The Hamiltonian 
should be divided into the ground state Hamiltonian and the residual interaction (and, 
if necessary, the external field). We consider without p-odd terms, hence = and as 
a consequence /° = 0. It is natural to take A° real, i.e. A° = 0, k° = and A^ = 2A°, 
Linearizing ([6]), we arrive at 



tnjJd{p,r)W6U = |rf(p,r) 



ih 



{{W, hl}5f^ + {W, hl}fl) + W{A15k_ - < (5A_ 



ih 



dt 



Jd{p,r)W6f. = Jd{p,r) 



-^{{{W,Al}}6^_-{{W,^l}}6A_ 
ih 



-{{W,hl}6U + {W,hl}f' 



ih , 



--i{W,Al}5K^-{W,Kl}5A^) 



ih^Jd{p,r)W5K+ = Jd{p,r) 



whl6K^ + y {ly, hijnl - W6 A_{i - /°; 



ih 



dt 



A° - ^{{{W, hl}}6K^ + {{W, /° }}5A_) 
d{p, r)WSK_ = J d{p, r) [whl6K+ + Whln^ - W6 A^l - /°) + WAlSf^ 
-^{{{W, hl}}6K^ + {{W, hl}}Kl + {{W, Al}}5U + {{W, fl}}6A^ 



(12) 



Let us look more closely at the variation of the gap, (5A(r, p), which should be expressed in 
terms of 6K,{r,p). In the framework of the method of moments, this can be done quite easily. 
According to (Q, the relevant variations are connected by the integral relation 

d^p' 



<5A(r,p) 



-v{\p-p'\)6K{r,p'). 



(13) 



(27r^)3 

Substituting this into ( JT2l) and changing the variables p ^ p', we obtain immediately the 
desired result. Examples for this type of calculation can be found in Appendix A. 

Until this point, our formulation is completely general. To proceed further, we are forced to 
make some approximations to get rid of higher-rank moments and to obtain a closed set of dy- 



namical equations for the second-rank moments. The usual problems of the method of moments 
are connected with integrals of the type Jd{p,r)WA^SK, /(i(p, r){iy, A°}5k, /(i(p, r)H^/°(5A, 
/ d{p, r){W, A^}6f etc. The variations 6K,{r, p), 6f{r, p) are integrated here with weights which 
are more complicated functions of r, p than the simple functions W. So the problem arises: 
how to express these integrals via the moments / d{p, r)W6K{r, p), / d{p, r)W6f{r, p) which we 
work with? To solve this problem, one can develop the functions /''(r,p), K°(r,p), A°(r, p) in 
a Taylor series around some point (f , p). However, another problem appears on this way: what 
to do with the higher-order moments which will inevitably be generated by the Taylor series? 
We suggest to neglect them, because it is natural to expect that the influence of higher-order 
moments on the dynamics of lower-order moments will be small [7] . 

A few words about the choice of the point (f , p). It is known that all the dynamics happens 
in the vicinity of the Fermi surface. Therefore, the choice of the momentum p is obvious: 
it should be equal to the Fermi momentum PF(f). The choice of f is more complicated, 
because it depends on the nature of the mode under consideration. For example, in the case 
of surface vibrations, f should be taken somewhere near the nuclear surface R. In the case of 
compressional modes, it is more appropriate to take r somewhere inside the nucleus. In any 
case, it is a rather delicate problem and every particular example requires the careful analysis. 
In principle, the value of r can be used as a fitting parameter. 



3 Equations of motion 

The model considered here is a harmonic oscillator with quadrupole-quadrupole residual inter- 
action [T]. The corresponding mean- field Hamiltonians for protons and neutrons are 

h^ = ^ + \^^"^" + j:y^TBl,{t)q,,{v) - /.^ (14) 

where B^^ = KnnQ% + i^npQl^, Bl^ = KppQl^ + HnpQl^, q2p,{j) = ^167t/5r^Y2^{9,(f)), /x^ being 
the chemical potential of protons (t=p) or neutrons (r=n). Krr' is the strength constant, and 
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we suppose Knn = i^pp- Qln is a component of the quadrupole moment 

Ql.it) = I d{p, r)g2^(r)r(r, P,t) = J dip, r)q2,iv)[f5ir, p) + ^^(r, p, t)] = + 5Ql^{t). 

(15) 

The Hamiltonian is divided into the equihbrium part and the variation h\\ 

= E {-ir5Bl^{t)q,,{v), (16) 

where we took into account that in an axially symmetric nucleus Q^ii ~ ^$±2 ~ 0- Obviously, 
in this model /iL = and hZ^ = 2h'^. 

The required Poisson brackets together with the definition of and cjj are written out 
in appendix B. For the sake of simplicity, we will neglect in (IT^ all terms proportional to h'^ 
(quantum corrections) except the simplest one, namely the term with h^, which we will keep 
in order to have an idea about the possible infiuence of quantum corrections. We assume that 
A° does not depend on the direction of p, i.e. A''(r, p) = A(r,p) (see Appendix A). Let us 
introduce the collective variables 

Q±(t) = J c/(p, r)xz5fl{r, p, t), Ql{t) = J c/(p, r)xz5Kl{r, p, t), 

P±{t) = Jd{P:r)p^Pz5flir,p,t), P±(t) = Jd{p,r)p^pjKl{r,p,t), 

Llit) = I dip, v)L5fl{v, p, t), Llit) = J dip, r)M<(r, p, t), 

Il{t) = jd{p,v)iy5fl{v,p,t), il{t) = jd{p,v)iy5Kl{v,p,t). (17) 

By definition, = = P_ = Q_ = 0. The analysis of the dynamical equations shows that 
L+ = /_!_ = I/_ = /_ = 0, and we are left with the following eight equations: 

1. ihQ^ = ih — L_, 

m 

2. zhQ^ = {2ho-\Vo\I^^)Q^ + '-^^-^Ia,L., 

3. ihQ_ = (2/io- |K,|//f)g+ + 2A;4Z(t) +2A°Q+, 



5. zhP+ = {2ho-\Vo\ll^)P^-'-^\Vo\l9ML^-u;U- 

8 



^2 

6. ihP. = (2/io-|K)|//p^)n + 2A°P+-yA;oZ(t), 



m n 



7. ihL_ = ih 

m 

8. thi^ = -in[mujlQ+ + 2{{z^) - {x'^))Z{t)]. 



Here Z^{t) = 12^1^(1) (see Appendix B) with pf^ = 1[kQI + RQl] and (3f^ = ^[kQI + 
kQI], (x^) = /c/(p,r)xVo(r,p), {z^) = /c/(p, r)zVo(r, p), h = 2 /rf(p, r)K°(r, p), h = 
2jd{p,r)x^z^K°{r,p) and ul = lu^ ± ul The integrals J/i^, J/p^, J^^f, Ji;^>, hp and 
are defined in Appendix A. All coefficients should be calculated at the point (f, pi?(f)), for 
example, A° = A°(r, p^(f)). Note that the first equation does not contain A-dependent terms. 
It has the typical structure [U [7] which is characteristic for coordinate moments of the continu- 
ity equation. The last equation does not depend on pairing parameters, either. However, this 
is due to the symmetry property of the operator / (see fl77l) ) and has nothing to do with the 
continuity equation. 

Let us compare these equations with the respective equations of [3], which were derived 
using the approximation 

A° = const, 6A{r, p, t) = 0, (19) 

resulting in a violation of the continuity equation. Written in terms of the variables (fT7|) . they 
read 

ihQ+ = i^— L_ + 2A°g_, 
m 

ihQ_ = 2A°Q+ + 2kiZ{t), 
m 
~2 

ihP. = 2A°P+ - —koZ{t), 

ihL. = ih\—P+-mujlQ+-2((z^) + (x'^))Z(t)], 
m 

ini_ = -in[mujlQ+ + 2{{z^) - {x^))Z{t)]. (20) 

The most evident difference between the sets of equations f|T8|) and fl20l) is the number of 
equations: eight and six, respectively. How could this happen? The approximation f|T9|) . used 
in [3], makes the third equation of f|T2|) trivial (i.e. its right-hand side becomes equal to zero 



ihP, = -ih—lculL^ - Lj^IJ + 2A°P_, 



identically), producing two integrals of motion = and P+, which should be included in 
the set (PU]). 

The second and the most important difference between the two systems of equations con- 
cerns the first equations of (fT8l) and (l20l) . The dynamical equation for the variable Q+ in (120|) 
contains the additional (in comparison with (flSll ) term 2A°Q_, whose existence is a direct 
consequence of the violation of the continuity equation. This is the only place where the vio- 
lation of the continuity equation appears explicitly. There are more differences between f fTSj) 
and (l20l) . which are all connected with the approximation f|T9|) . For example, the dynamical 
equation for in f[T5]) has the term \Vo\Ip^P- instead of the corresponding term 2A°P_ in 
(I20p . The integral 1^^ contains the contributions from 6k (the factor 2A°) and from 6 A as well 
(see Appendix A, formulae (H2|) and (!57|) ). 

And the last difference: the set of equations (120!) has the pleasant property that its eigen- 
modes can be found analytically, contrary to those of the set f llSp . There is no necessity to 
explain how important and convenient it is to have (even approximate) analytical solutions of 
a problem. It turns out that one can find an approximation which allows one to get analytical 
solutions of ffTS]) without violating the continuity equation. 

From general arguments one can expect that the phase of A (and of k, since both are 
linked according to equation (Q) is much more relevant than its magnitude, since the former 
determines the superfluid velocity. After linearization, the phase of A (k) is expressed by 6A_ 
{Sk_), while 5A+ (5k+) describes oscillations of the magnitude of A (k). Let us therefore 
assume that 

(5fi:+(r,p) < 5K_(r,p). (21) 

This assumption was explicitly confirmed in ref. [S] for the case of superfluid trapped fermionic 
atoms, where it was shown that 5A_|_ is suppressed with respect to 6A_ by one order of A/Ep, 
where Ep denotes the Fermi energy. 

The assumption fl2T|) does not contradict the equations of motion and allows one to neglect 
all terms containing the variables (5+ and -P+ in the equations No. 3, 6, and 7 of flTSl) . In this 
case the "small" variables Q+, P+ will not affect the dynamics of the six "big" variables Q+, 
P+, L_, /_, Q-, P-. This means that the dynamical equations for the "big" variables can be 
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considered independently from that of the " small" variables, and we will finally deal with a set 
of only six equations. Adding the isospin index r =(n,p), we have 

i?iQ\ = ih — L^, 
m 

ihQ_ = 2A°^Ql + 2klZ^{t), 
ihP_ = 2A°^ PI - —k^z^it), 

iht = -ih[muj:_^Q^^ + 2{{zy - {xY)Z^{t)]. (22) 

It is interesting to compare the results of the two different approximations: the sets of equations 
(122|) and (120|) . The difference is minor: the factor |Vo|/pjf in the third equation of fl22|) instead 
of 2A° in ( |20|) and the absence of the term 2A^Q_ in the first equation of (122|) contrary to (!20!) . 
Calculations show that numerically the factor |Vo|/pjf^ is not so far from 2A° (see fig. 1), so one 
can conclude that the approximations fl2Tl) and f|T9l) lead to similar dynamical equations. On 
the other hand, the approximation fl2T]) is undoubtedly better than f|T9|) . because it is physically 
motivated and it does not violate the continuity equation. 

In this paper we will treat the reduced set of equations ( 122|) that can be solved analytically. 
This allows us to compare our results with those of [3J and to assess the quantitative effect of 
the correct treatment of the continuity equation. The numerical analysis of the full problem 
(lisp with all coefficients calculated in the microscopic approach will be postponed to a future 
publication. 

4 Analytical solution 

First of all, we rewrite the equations fl22|) in terms of the isovector and isoscalar variables 
Q^ = Ql~ Ql, Q± = QI + Ql, Q^ = Ql- Ql, Q± = QI + Ql, and so on. In order to 
separate the isovector and isoscalar sets of equations, we employ the standard approximation 
which works very well in the case of collective motion [7|: 

Ql/N = ±Ql/Z, (23) 
11 



where {Z) is the number of neutrons (protons) and the sign + (— ) is utihzed for the isoscalar 
(isovector) motion. 



4.1 Isovector eigenfrequencies 

The set of equations describing isovector excitations reads 





= ih — L , 
m 


iTiQ 






3 


ihP_ 


= 2A'P_^--xoKQ+: 


ihL_ 


2 2 

= ih — P , — iTi2muj ijjQ , 
m — + 


ihl_ 


= -a)(l -02mcu 



(24) 

The following notations are introduced here: Xo = 12ko, ip = [{I + |)(1 — a — ^) — 
k'. = akf + (z/ - Tx)kr^ kf = kf ± fcf, i = 0,4, A' = z/A" + ttA?, A' = z/A'^ + ttAp, 



2A- = m{i;^)\ V = N/A, vr = Z/A, A = N+Z, ^ = {u-7r){Q-,-Ql,)/Qoo, Qoo = QSo+Qoo- 
We use the standard [lOj definition of the deformation parameter 6 = 3Q2o/(4Qoo)5 where 
QoQ = Jd{p,r){x^ + + z^)f{r,p) = 2{x^) + (z^) and Q20 = 2(z^) - 2{x'^) are monopole 
and quadrupole moments, respectively. Usually one takes the isovector strength constant ni = 
|('^nn — i^np) proportioual to the isoscalar one, kq = |(fi;„„ + Knp), i-e., = a/to, a being a 
fitting parameter. Following ref. [T], we take a = —2. For the isoscalar strength constant we 
take the self consistent value (see appendix B). 

This set of equations has two integrals of motion: 



mA' 

and 



^ - 2mu^i: 
4mA' 



ih 



XoK + 2A 



-Q — L_ = const (25) 



1 — £ 

L_ + ih2muj'^5{l - a) — - — ^-rQ = const. (26) 

Xok'4^ + 2A'— 



By definition, the variable Q„ is purely imaginary because is the imaginary part of the 
anomalous density k. Therefore Eq. fl26l) implies that the relative angular momentum /_ 
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oscillates in phase with the relative quadrupole moment Q_ of the imaginary part of k. Anal- 
ogously, one can interpret Eq. (!25|) saying that the variable L_ oscillates out of phase with the 
linear combination of two variables Q_ and P__ describing the quadrupole deformation of the 
anomalous density k in coordinate and momentum spaces, respectively. 

Imposing the time evolution via e*^*/^ for all variables, one transforms the equations 
into a set of algebraic equations, whose determinant gives the eigenfrequencies of the system. 
We have 

E^E^ - 2V^E^ + SA'A'nW^ + 4{hoo)^S\l - a){l - 0(1 - «0 - Xo^'Kh^/m'} = 0, (27) 

where = 2 A' A' + h'^iu^[{l + |) (2 - a - - pa^]. The solution E = corresponds to the 
integrals of motion (125!) and (126|) . Two nontrivial solutions of (1271) 

El=V^± ^Vl - SA'A'h^io^ilj - A{hoo)'^6^{l - a){l - 0(1 - O + Xo^'K^V^'^ (28) 

describe the energy E^ of the Iso Vector Giant Quadrupole Resonance (IVGQR) and the energy 
E_ of the scissors mode. It is worth noting that contrary to the case without pairing [I| the 
energy E^ does not go to zero for deformation 6 = 0. However, this does not contradict 
the known quantum mechanical statement that the rotation of spherical nuclei is impossible.lt 
is easy to see from fl26l) that the relative angular momentum /_ is conserved in this case, 
/_ = const, so this mode of a spherical nucleus has nothing in common with a vibration of 
angular momentum. The calculation of transition probabilities (see below) shows that it can 
be excited by an electric field and it is not excited by a magnetic field. Our estimate of the 
energy of this mode gives a value of about 2.88 MeV, which is not far from the result of M. 
Matsuo et al [TT], who studied the isovector quadrupole response of ^^^Sn in the framework of 
QRPA with Skyrme forces and found the proper resonance at ~2.2 MeV. 

It is known [H |2] that without pairing the scissors mode lies at non-zero energy only due to 
Fermi Surface Deformation (FSD). Let us investigate the role of FSD in the case with pairing. 
Omitting in (12^ the variable P_,_ responsible for FSD and its dynamical equation, we obtain 
the characteristic equation 

E'{E' - 2(M'[(1 + ^)(1 - « - - = 0, (29) 
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which coincides with the analogous equation of [T] derived without pairing for ^ = (A^ = Z). 
The two solutions Ef^^ = and El^^f^ = 2{huj)'^[{l + |)(1 — a — ^) — |a^] demonstrate in an 
obvious way that the role of FSD is not very important for IVGQR, whereas it is crucial for 
the scissors mode and the ISGQR, whose energy in the approximation ^ = can be obtained 
from the IVGQR by assuming a = 1 (see below). 

4.2 Transition probabilities 

The transition probabilities are calculated with the help of linear response theory. The detailed 
description of its use within the framework of the WFM method can be found in [1] , so we only 
present the final results. 

Electric quadrupole excitations are described by the operator 

F = Fi^ = j2kis), h, = er^Y,,. (30) 

The transition probabilities of the isovector modes are 

Bm. ^ 21 < .l^a 10 > P ^ QS„ ^ '1'^^^ - ^^f > - f - . (31, 

Magnetic dipole excitations are described by the operator 

^ . . eh 

^ = ^im = E/im(^), /iM = -^V(ryi^) ■ [r X VW, f^N = (32) 

Their transition probabilities are 

B{M1)^ = ^—{l-a){l-OQoo^ E [E-^-V ] ^ ^ 

Multiplying the B(M1) factors of both states by their respective energies and summing up, 
we find the following formula for the energy-weighted sum rule 

— 2 

EscB{Ml)sc + E,,B{Ml)i, = (1 _ «)(i _ ^) QP^s^^l, (34) 

ZTT 

The same manipulations with the B(E2) factors give 

5 

EisLBiE2)jsL + E,sB{E2),s + E,,B{E2),, + E,,B{E2)i, = e''-—QU^ + S/S). (35) 

m /TT 
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It turns out that both expressions coincide exactly with the respective sum rules calculated in [T] 
without pairing. Does this mean that there is no contribution to the sum rules which comes from 
pairing? Of course not, because the value of the mean square radius Qoo should be calculated 
with the ground state wave function which depends on pair correlations: Qoo = Z]j '^^^(^k^K)• 
This is a good place for discussing the deformation dependence of the energies and transition 
probabilities of the isovector modes. First we recall the relevant formulae without pairing [1]: 




B(Ml,;^JL».W^;|^^|y^.^. (36) 

where 63 = 1 + 6/3, the superscript "0" means the absence of pairing and we assumed a = —2. 
For the sake of simplicity we put ^ = 0. The scissors-mode energy is proportional to 6, which 
becomes evident after expanding the square root: 



3_ / 3^6^ 9 S 
+16 51+128 53 

At a first glance, the transition probability, as given by formula (136|) . seems to have the 
desired (experimentally observed) quadratic deformation dependence. However, due to the 
linear 5-dependence of the factor Eg^ in the denominator, the resulting 5-dependence of B{M1)^^ 
turns out to be linear, too. The situation is completely different when pairing is included. In 
this case, the main contribution to the scissors mode energy comes from the pairing interaction. 
Esc is not proportional to 6, and the deformation dependence of B{Ml)sc becomes quadratic in 
excellent agreement with QRPA calculations and experimental data [2l [T2 | [T3 t [T ^ [T5l [T6l [T7t [T8] . 

The deformation dependence of B{Ml)iy is quadratic in 6, even without pairing, because 
the energy Eiy is not proportional to 6 and depends only weakly on it. The inclusion of pairing 
does not change this picture. 

4.3 Isoscalar modes 

The set of equations describing isoscalar excitations reads 

ihQ^ = ih — L_, 
m 
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i^Q- = {xo[kt + a{u-7r)k^]+2A'}Q+, 

ihP+ = 2A'P.-ihmiu^{[l + -{l-aO]L--6{l-aOl-}, 

3 

ihP. = 2/\'P+-—xo[kt + <^{y-TT)k^]Q+, 

2 2 
ihL_ = ih—P, + ih2miu^(l + -S)a^Q+, 
m 3 

ihi_ = 0. (38) 

As it is seen from the last equation, the angular momentum J_ = J° + P is conserved, as it 
should be, since we work with the rotationally invariant mean field Hamiltonian iHM . 
Assuming for simplicity N = Z, we find the following characteristic equation 

E^{E^ - E'^[4AA + 2e^] - xo^koh^/m^} = 0, (39) 
where = h^uj'^{l + |). The two solutions of ([39]) 



El = 2AA + e^± V(2AA + e^y + xoAkot/m'^ (40) 

give the energy E^ = Eis of the IsoScalar Giant Quadrupole Resonance (ISGQR) and the 
energy E^ = Ejsl of the IsoScalar Low-Lying Excitation (ISLLE). If one neglects the quantum 
correction (the term with ko), one finds for the ISGQR energy the expression 

E'gqr = 2e' + 4AA, (41) 

which is reduced to the standard value Eqqr = V^huj{l + |5) in the case A = 0. In this 
case, the energy of the low-lying mode disappears. The transition probabilities of the isoscalar 
modes in the approximation N = Z are obtained from (|3T1) by taking a = 1. 

It is important to note that the very existence of the ISLLE relies on two factors: 1) pair 
correlations and 2) quantum correction. With the parameters given in section 4.4, its energy 
and transition probability for ^^^Dy can be estimated to be Ejsl = 1-0 MeV, B{E2) = 41.4 
W.U.. These numbers are of the right order of magnitude [1]. Nevertheless, we do not dare 
to compare them with an experiment until all terms proportional to h"^ in f|T2|) are taken into 
account. The accurate calculation of the quantum correction and the comparison of Ejsl and 
B{E2)jsL with experimental data will be postponed to a future publication. 
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4.4 Numerical results for the scissors mode 

We have reproduced all experimentally observed qualitative features of the scissors mode. We 
understand that the harmonic oscillator model with QQ residual interaction is too simple to 
give a precise quantitative description of the experimental results. Moreover, even within this 
simple model we had to make the additional approximation fl2T|) . Nevertheless let us calculate 
the energies and B{M1) factors to get an idea of the order of magnitude of the discrepancy 
with experimental data. We will also compare our results with those of [3J in order to see the 
effect of the violation or non-violation of the continuity equation. 

Results for most of the nuclei where this mode has been observed are presented in Table [1] 
and in Figures 1-3. The formulae fl28p and 0331) were used with the following values of the 
parameters: a = -2, Qqq = A^R^ R = roA^/\ ro = 1.2 fm, uj^ = uj^/[{l + - |5)V3]^ 

hujo = 41/A^/3 MeV, h'^/m = 41.803 MeV fm^. The gap A as well as the integrals I^^ 
and kf) were calculated with the help of the semiclassical formulae for n{r,p) and A(r,p) (see 
Appendix A), a Gaussian being used for the pairing interaction with rp = 1.9 fm and Vq = 25 
MeV. The dependence of A(r,pi?(r)) and Ip^{r,pF{r)) (calculated for ^^"^Er) on the coordinate 
r is demonstrated on fig. 1. We checked that the dependence of A{r, pf{t)), Ip^{r,pF{r)) and 
Pi?(r) on the direction of the vector r is negligible. The best agreement of the theory with the 
experiment is obtained at the point f where the integral Ip^{r,pF{r)) has its maximum. The 
values of fc^, vary smoothly from kQ = 62.5 for A = 134 to /cq" = T2.2 for A = 196. The analysis 
of Table 1 shows that the overall agreement of the theoretical results with the experimental 
data is reasonable. It is, of course, not perfect, but the influence of pairing, especially on the 
B{M1) values, is impressive. Without pairing, the calculated energies (column A = 0) are 1 - 
2 MeV (1.5 - 2 times) smaller than Eexp, and the B{M1) factors (column A = 0) are 3-7 times 
larger than the experimental values. The inclusion of pairing changes the results (columns 
new) drastically: the discrepancy between calculated and experimental energies is reduced to 
5 - 20% and the calculated transitions probabilities are reduced by a factor of 1.5 - 3. The 
influence of fco is negligibly small, being of the order of ~ 2%. 

The influence of the correct treatment of the continuity equation can be estimated by 
comparing the columns new (results with the fulfilled continuity equation) and old (results 
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Table 1: Scissors mode energies Esc (in MeV) and transition probabilities B{Ml)sc (in units of 
exp: experimental values, old: old theory |3j with violated continuity equation, new: new theory with 
fulfilled continuity equation, A = 0: theory without pairing. The experimental values of Egc, S and 
B{M1) are from Ref. [19] and references therein. 
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Figure 1: The pair field A(r,pi?(r)), the value of A = 2A = \Vo\Ip^{r,pF{r)) and the Woods-Saxon 
potential Vwsif) the functions of radius r. 

from ref. [3] with the violated continuity equation). It turns out that the correct treatment of 
the continuity equation has a substantial effect on the results: on the one hand, it leads to a 
decrease of the scissors-mode energies by 0.8 - 0.9 MeV (in comparison with the old results), 
improving the agreement with the experimental data, but on the other hand, it results in an 
increase of the transition probabilities by a factor 1.4 - 1.5, deteriorating the agreement with 
experiment. 

What can be done to improve these results? An obvious idea is to get rid off the approxi- 
mation fl2T|) . i.e., to solve the full set of eight equations f|T8l) . calculating all integrals within a 
microscopic approach. The next possible step is to perform a self-consistent calculation with a 
more or less realistic interaction. 

Another point which should be clarified is the role of the spin-orbit interaction. It is known 
P] that the experimentally observed low-lying magnetic dipole strength consists of two sepa- 
rated parts: orbital excitations in the energy interval ~2 - 4 MeV and the spin-flip resonance 
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Figure 2: Scissors mode energies as a function of the mass number A for the nuclei listed in Table 1. 
Enew- new theory, Eoid- old theory [3], Eq: theory without pairing. 
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Figure 3: Scissors mode transition probabilities B(M1) as a function of the mass number A for the 
nuclei listed in Table 1. Bnew' new theory, Bold- old theory Bq: theory without pairing. 
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ranging from 5 to 10 MeV excitation energy. So, for the full description of the scissors mode it 
would be important to consider also the spin degrees of freedom. 

5 Conclusion 

In conclusion, we presented a generalization of the method of Wigner function moments which 
allowed us to include pair correlations without violating the continuity equation. The method 
was exemplified by the calculation of isovector and isoscalar giant resonances and low-lying 
excitations in the harmonic oscillator model with quadrupole-quadrupole residual interaction. 

The analytical formulae, derived in a slightly simplified model (approximation fl2T|) ). repro- 
duce very well the experimentally observed deformation dependence of the energy and B{M1) 
factor of the scissors mode. We performed calculations for most of the nuclei where this mode 
has been observed, and our results are in reasonable quantitative agreement with the experi- 
mental data, the pair correlations being extremely important. 

The joint action of pairing and quantum corrections lead to the appearance of a low-lying 
isoscalar excitation. The low-lying isoscalar, as well as isovector, modes are very sensitive to 
the parameters of the interaction and to all details of the ground-state input. Their study on 
the basis of the full (8 variables) set of dynamical equations f|T8l) . with normal and anomalous 
densities calculated within a microscopic approach, are in progress. 
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Appendix A 

Let us consider integrals containing 5k, and SA. 

The relevant term of the first equation in f|T2|) with W = PxPz reads 




(42) 
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Substituting A and SA according to the expressions and f[T^ . one gets 
Changing in the second part of the integral the variables p ^ p', one finds 
Let us analyze in detail the integral over p' 



Now it is necessary to clarify the p dependence of hP . There are good reasons to assume that 
it depends only on |p| = p, i.e., it does not depend on the angles in the momentum space. In 
fact, according to the semi classical formula [5] 

.(r,p) = l , ^(^'P^ (46) 
V/^'(r,p) + A2(r,p) 



the anomalous density can depend on the direction of p only via A(r, p). However, as can be 
seen from the semiclassical gap equation [5] 

A(r,p) = — /- — TT^v(\p-p\) , (47) 



(27r/i)3 y^/i2(r,p') + A2(r, p') 

there are no reasons to introduce the angular dependence (in momentum space) of A as long 
as the Fermi surface is spherical. So, it is quite natural to take A(r,p) = A(r,p) and, as 
a consequence, ^(r, p) = k{v,p). In this case the integral /p, can be integrated over angles 
analytically. To this end we expand the pairing force in a standard way [2TI [5]: 

oo 

v{\p - P'l) = E <P^P) E Vm{0. mUd', 0')- (48) 

1=0 "I 

Having in mind that p^Pz = \/|f^P^(^2-i — ^21) we find 

2p.P.J^^,>^\r,p')vo{p,p') - 2^[Y,_,{9,^) - Y,,{9,^)] J ^^,'^'i^^P')Mp,p'). 

(49) 
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If we use a Gaussian as pairing interaction [S]: 

^^(|p-p1)=/5e-"l^-P'l% (50) 
with (3 = —\Vo\{rp^/7^Y and a = r'^/Afi^^ the expansion coefficients of the force read [21] : 

vi{p,p') = An/Se-'^^f'-^P'VU-ix), (51) 
where ji{—ix) is the spherical Bessel function and x = 2app'. We need the ffist three coefficients 
vo{p,p') = 47r/3e-"(^'+^"Vo(x), 0o(x) = -sh{x) = {1 + ^ + ...), (52) 

v,{p,p') = 47r/?e-"(^'+^")a;0i(x), 0i(x) = ^[ch{x) - -sh{x)] = 1(1 + ^ + . . .), (53) 

X 3 10 

V2{p,p') = 47r/3e-"(p'+^'")xV2(x), 02(x) = ^[(1 + ^)sh{x) - -ch{x)] = 1(1 + ^ + . . .). 

(54) 

With the help of these expressions one finds 



p,p,|VolC(r,p), (55) 



pp 

op 

where 

3 

a^^P) = / '^°(^'^'') h(2app') - 4aVV2(2app')] e-^^'^'c^p'- (56) 

Substituting ( l55l) into (jSj) one gets finally 

5 = |m|rf(p,r)j;p^(r,p)p.p,5/t_(r,p,t) 

^ \Vo\I;^{v,pf) Jdip,T)p,pjK.ir,p,t) = \Vo\i;^if,PF)P-it)- (57) 
The third and fourth equations of ( fT2l) contain the terms 

Tf = |rf(p,r)iy(l-/°)5A±. (58) 

As above, we assume that A does not depend on the direction of p. On the basis of the 
semiclassical formula 

l-/°(r,p)= ^ (59) 
' //.2(r,p) + A2(r,p) 
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we can assume /^(r,p) = /^(r, p), in agreement with the assumption of a spherical Fermi 
surface. Calculating (!58|) with W = xz, PxPz, L and ly we obtain 

= (r,p^)Q±(t), = \Vo\ll^{v,PF)P±{t), (60) 

Tt = ml^{v,pp)L^{t), Ti = |Fo|/4^(f,PF)/±(t), (61) 

where 

3 

//f (r,p) = ^^^'""'Jl^ - fl{r,p')]M'2app')e-^^"p''dp', (62) 

= i;^^""'' /[I - /^(r,p')]0i(2«PP')e'"^'V''rfp'- (63) 

7 

(r,p) = ^^^""'Vf^ " /°(r,p')]02(2appOe-"^'W, (64) 
The second and third equations of f|T2l) contain the terms 

= 2 j ci(p, r) ({VT, A°}(5fi:+ - {IV, /t°}(5A+) (65) 

and 

= /rf(p,r){H^,A°}<5/_, (66) 
which require the knowledge of derivatives of A°(r, p). They are found with the help of formula 



The derivative with respect to p reads: 

9A0 IK) I 



Idp{r,p)pi, (67) 



dp^ 2h' 
where 

l9piT,p) = -^e-°^' / K%v,p') k(2app') - 2ap'^M'2app')] e^^^^p^dp'. (68) 

To calculate the derivative with respect to r, we approximate /t°(r,p) by formula ( l46l) . As a 
result, we obtain an integral equation for dA^{r, p)/dxi, with a kernel which is strongly peaked 
at p' = pf, which allows us to simplify the equation by replacing dA'^{r, p')/dxi under the 
integral by dA^{T,pp)/dxi. Finally we have 



gA°(r,p) 

dxi 



f3{r,p,pF)muJiXi, (69) 
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where 
and 



{2'KhY 



Mp^p') - Mpf,p') 



1 + 7(r,PF) 



h{r,p')A{r,p') 



7(r,p) 



:Mp^p') 



h'^{r,p') 



{2TiTif ' 2[/i2(r, p') + A2(r, p')fl'^ ' 

Calculating now fl66l) with W = xz, we find 

r BA^ BA^ \V I 

G^^ = 2j d{p, r) + x^)5f.{v, p, t) = -L^I,^(r,pp)L_{t). 

The calculation of (!66i) with W = PxPz gives 



where 

Idr{r,p) 



ox oz 



)(5/_(r,p,t) = \V^\hr{f,PF)HL^{t)-ujtl-{t)l 



mr„ 2 
__E_g-"P 



)(2app) 



/;,(r,p')A(r,p') - h'^{r,p')p{T,p',pf) 



[/i2(r,p') + A'(r,P')]^/' 



Calculating (!65|) with the weight xpz, one gets 
with 



1 /€°(r,p')[0o(a:) - Aap'^M^) + 4:ay%{x)]e-''P"p'^dp\ 



where x = 2app'. It is evident that K^^'' = K^^"^ . As a result 

i^^ = -2^Jif>(f,p^)P+(t), K' = Q. 



Appendix B 

The necessary Poisson brackets are 
{xz, hi) = ^L, {pxp„ hi} = -m{ul^zp,x + ujl^xp^) = -^[(^f + ujI^)L - {ul^ - 
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{L, hi} = —PxPz - m{ujl + ul )xz, {ly, h^} = m{ujl - luI )xz, {xz, hi} = 0, (79) 
{PxPz, hi} = 26Bl,{t){xp, - 2zp^) - Q[(3l,{t) - f3^,,{t)]xp, 

-12[Pl,{t)izp, + xp,) + PUt)yp, + P;sit)yP^l (80) 

{L, hi} = -25Bl^{t)xz - l2[l3Ut){z^ + x^) + I3l^{t)yz + PUt)yx] - (t) - I3l^{t)]xz, (81) 

{4, hi} = Q5Bl,{t)xz + l2[[5lS)yx - PUt)yz - l3lS){z' - x^)] - 6[/3[i(t) - ^Ut)]xz, (82) 

{{W. K}} = 0, {{^z, hi}} = 0, hi}} = 24/3[3(t), {{L, hi}} = 0, {{4, hi}} = 0. (83) 
Here 

^^=KJl^ + Rjf^, P^ = nJ^ + Rj-, (84) 

J^.(t) = |d(p,r)x.Vr(r,p,t), J[3W = ^Q;W, 5i?2^oW=2/5l3W-/5nW-/522W- (85) 

Using the definition [10\ Q20 = ^Qoo^ and the formula ^20 = — a;^ — y^, one finds from f|T6|) 



^^' = ^'(1- A^2o), a;^' = a;'(l-^i?2o)- (86) 



These oscillator frequencies depend on the strength constants and Knp. In the isoscalar case 
the constant kq = |(/tnn + /tnp) is fixed by the self-consistency condition [10] uj^x"^) = ujy{y^) = 
uj1{z'^)^ which allows one to find for hq and the oscillator frequencies of the isoscalar field the 
following expressions |20j : 

«:o = -^, u:l = \{uf + uf)=u:\l + U\ = \{uf + u:f) = u:\l - (87) 

where = cuV(l + U) = + |5)i(l - |(5)i], and huo = Al/A^/^MeV. 
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